library(Signac)
library(Seurat)
library(tidyseurat)
options("restore_Seurat_show" = TRUE)
library(tidyverse)
library(ggpubr)
source('00_util_scripts/mod_bplot.R')

# import peak mtx ----------
atac.norm <- Read10X_h5('~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/filtered_peak_bc_matrix.h5')
atac.norm[1:5,1:5]

atac.norm |> glimpse()

atac.meta <- read_csv('~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/singlecell.csv') |>
  column_to_rownames('barcode')

chrom_assay <- CreateChromatinAssay(
  counts = atac.norm,
  fragments = '~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/fragments.tsv.gz',
  sep = c(":", "-"),
  min.cells = 10,
  min.features = 200
)

aobj <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = atac.meta
)

aobj[['peaks']]

granges(aobj)

atac.psor <-
  Read10X_h5('~/append-ssd/alaria2/sun2024psoria/psoria_fullq/outs/filtered_peak_bc_matrix.h5')

atac.meta <- read_csv('~/append-ssd/alaria2/sun2024psoria/psoria_fullq/outs/singlecell.csv') |>
  column_to_rownames('barcode')

chrom_assay <- CreateChromatinAssay(
  counts = atac.psor,
  fragments = '~/append-ssd/alaria2/sun2024psoria/psoria_fullq/outs/fragments.tsv.gz',
  sep = c(":", "-"),
  min.cells = 10,
  min.features = 200
)

aobj.psor <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = atac.meta
)

# keep only standard chromsomes
peaks.keep <- seqnames(granges(aobj.psor)) %in% standardChromosomes(granges(aobj.psor))
aobj.psor <- aobj.psor[as.vector(peaks.keep), ]

# add genome annotation --------
## extract gene annotations from EnsDb packages
## create from gtf 
# annotations <- GetGRangesFromEnsDb(ensdb = ensdb.v112)

library(AnnotationHub)
ah <- AnnotationHub()

# use latest ensembl
latest.ensdb <- query(ah, "EnsDb.Hsapiens") |>
  tail(1)

# or the version same as cellranger-atac index
# latest.ensdb <- query(ah, "EnsDb.Hsapiens.v98")

# download ensdb may cost ~10 min
# load from cache cost 40s
annotations <- GetGRangesFromEnsDb(ensdb = latest.ensdb[[1]])

seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg38"

# add the gene information to the object
Annotation(aobj) <- annotations

Annotation(aobj.psor) <- Annotation(aobj)

# QC ------
# compute nucleosome signal score per cell
# compute TSS enrichment score per cell
# cost 1m40s
aobj.psor <- aobj.psor |>
  NucleosomeSignal() |>
  TSSEnrichment()
  
# compute blacklist frac via Signac 
aobj.psor$blacklist_fraction <- aobj.psor |>
  FractionCountsInRegion(blacklist_hg38_unified)

# add blacklist ratio and fraction of reads in peaks
aobj.psor <- aobj.psor |>
  mutate(pct_reads_in_peaks = peak_region_fragments / passed_filters * 100,
         blacklist_ratio = blacklist_region_fragments / peak_region_fragments)

aobj.psor |>
  DensityScatter(x = 'nCount_peaks', y = 'TSS.enrichment',
                 log_x = TRUE, quantiles = TRUE)

aobj$high.tss <- ifelse(aobj$TSS.enrichment > 3, 'High', 'Low')
TSSPlot(aobj, group.by = 'high.tss') + NoLegend()

aobj$nucleosome_group <- ifelse(aobj$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = aobj, group.by = 'nucleosome_group') +
  ggtitle('Nucleosome signal QC')

# qc violin plot
aobj.psor |>
  VlnPlot(
  features = c('nCount_peaks', 'TSS.enrichment', 'blacklist_fraction',
               'nucleosome_signal', 'pct_reads_in_peaks'),
  pt.size = 0
)

# remove cells of bad quality
# adjust threshold accordingly
# 5% < nCount_peaks < 95% * 2
aobj <- aobj |>
  dplyr::filter(nCount_peaks > 1000 &
           nCount_peaks < 30000 &
           pct_reads_in_peaks > 15 &
           blacklist_ratio < 0.05 &
           nucleosome_signal < 4 &
           TSS.enrichment > 3)

aobj.psor <- aobj.psor |>
  dplyr::filter(nCount_peaks > 1000 &
                  nCount_peaks < 4100 &
                  pct_reads_in_peaks > 15 &
                  blacklist_ratio < 0.05 &
                  nucleosome_signal < 4 &
                  TSS.enrichment > 3)

aobj

# Normalization & dim reduc & clustering -------
aobj.psor <- aobj.psor |>
  RunTFIDF() |>
  FindTopFeatures() |>
  RunSVD()

# LSI (TFIDF+SVD) component 1 will strongly correlate with seq depth
# So we will exclude it from downstream analysis
DepthCor(aobj.psor)

# SLM (smart local moving) clustering
aobj.psor <- aobj.psor |>
  RunUMAP(reduction = 'lsi', dims = 2:30) |>
  FindNeighbors(reduction = 'lsi', dims = 2:30) |>
  FindClusters(algorithm = 3)

aobj.psor |>
  DimPlot(label = TRUE, label.box = T, cols = DiscretePalette(36))

# gene activity matrix ---------
# cost ~5 min
gene.activities <- GeneActivity(aobj)

# add the gene activity matrix to the Seurat object as a new assay and normalize it
aobj[['RNA']] <- CreateAssayObject(counts = gene.activities)
aobj <- aobj |> NormalizeData(
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(aobj$nCount_RNA)
)

DefaultAssay(aobj) <- 'RNA'

FeaturePlot(
  object = aobj,
  features = c('KRT1','TRPV3','KRT14','CD8A'),
  pt.size = 0.1,
  max.cutoff = 'q95'
)

skin.marker <- list(
  kera = c('KRT1','KRT10','KRTDAP','KRT14','KRT5','COL17A1'),
  fibroblast = c('FBN1','LUM','COL3A1'),
  macrophage = c('CD163','MRC1','MARCO'),
  T.cell = c('CD3E','CD8A','IL17A','FOXP3'),
  DC = c('CD1A','CD207')
)

aobj |>
  DotPlot(c(skin.marker,'TRPV3'), cluster.idents = T, assay = 'RNA') +
  RotatedAxis()

all.marker <- aobj |>
  FindAllMarkers(assay = 'RNA',features = c(list_c(skin.marker),'TRPV3'),
                 only.pos = T) |>
  as_tibble()

all.marker |>
  dplyr::filter(p_val_adj < .05)

aobj <- aobj |>
  mutate(manual.main = case_when(seurat_clusters == 9 ~ 'T cell',
                                 seurat_clusters %in% c(1,11) ~ 'Macrophage',
                                 seurat_clusters %in% c(6,7) ~ 'Smooth muscle cell',
                                 seurat_clusters %in% c(2) ~ 'V3-lo KC',
                                 seurat_clusters == 10 ~ 'V3-hi KC',
                                 seurat_clusters %in% c(3,8,12) ~ 'Endothelial cell',
                                 .default = 'Fibroblast'))

aobj |>
  DimPlot(group.by = 'manual.main', cols = 'Paired') +
  ggtitle('scATAC-seq of healthy human skin')

publish_pdf('scatac.pdf', width = 140, height = 100)

## keep only fibroblast & KC -------
aobj.fk <- aobj |>
  dplyr::filter(manual.main == 'KC')

DefaultAssay(aobj.fk) <- 'peaks'

aobj.fk <- aobj.fk |>
  RunTFIDF() |>
  FindTopFeatures() |>
  RunSVD()

# SLM (smart local moving) clustering
aobj.fk <- aobj.fk |>
  RunUMAP(reduction = 'lsi', dims = 2:30) |>
  FindNeighbors(reduction = 'lsi', dims = 2:30) |>
  FindClusters(algorithm = 3)

aobj.fk |>
  DimPlot(label = T, label.box = T, cols = 'Paired')

aobj.fk |>
  FindAllMarkers(features = 'TRPV3', only.pos = T, assay = 'RNA') |>
  dplyr::filter(p_val_adj < .1)

fk.all.marker <- aobj.fk |>
  FindAllMarkers(only.pos = T, assay = 'RNA') |>
  dplyr::filter(p_val_adj < .05) |>
  as_tibble()

fk.all.marker |>
  group_by(cluster) |>
  slice_min(p_val_adj)

# Find differentially accessible peaks ------
# change back to working with peaks instead of gene activities
DefaultAssay(aobj) <- 'peaks'

## simple, acceptable & quick
da_peaks.wilcox <- aobj |>
  FindAllMarkers(only.pos = T, latent.vars = 'nCount_peaks')

head(da_peaks.wilcox)

## use this param will cost >2h for one cluster
da_peaks <- aobj |>
  FindAllMarkers(test.use = 'LR', latent.vars = 'nCount_peaks',
                 only.pos = T)

head(da_peaks)

## identify closest gene from granges
marker.peak.cluster <- da_peaks.wilcox |>
  as_tibble() |>
  dplyr::filter(p_val_adj < .05) |>
  pull(gene)

marker.peak.close0 <- aobj |> ClosestFeature(marker.peak.cluster)

region.v3 <- marker.peak.close0 |>
  as_tibble() |>
  dplyr::filter(gene_name == 'TRPV3') |>
  pull(query_region)

da_peaks.wilcox |>
  as_tibble() |>
  dplyr::filter(gene == region.v3)

der.v3hvl <- aobj |>
  FindMarkers(ident.1 = 10, ident.2 = 2)

deg.v3hvl <- der.v3hvl |>
  as_tibble(rownames = 'gene') |>
  dplyr::filter(p_val_adj < .05) |>
  pull(gene) |>
  ClosestFeature(object = aobj, regions = _)

deg.v3hvl |>
  as_tibble() |>
  filter(gene_biotype == 'protein_coding') |>
  filter(gene_name %in% c(kegg_mva, late.kc, 'TRPV3'))

# plot genomic region -------
pbmc |> CoveragePlot(
  region = region.v3,
  extend.upstream = 0,
  extend.downstream = 80000,
  assay = 'peaks'
)

# rds checkpoint ----------
aobj |> write_rds('mission/FPP/psoriasis/sun24.ctrl.signac.rds')

aobj <- read_rds('mission/FPP/psoriasis/sun24.ctrl.signac.rds')

# transfer anchor from scRNA-seq ----------
sobj_epi <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

transfer.anchors <- FindTransferAnchors(reference = sobj_epi, query = aobj,
                                        reduction = 'cca',query.assay = 'RNA')

predicted.labels <- transfer.anchors |>
  TransferData(refdata = sobj_epi$bill_fine,
               weight.reduction = aobj[['lsi']],
               dims = 2:30)

predicted.labels|> head()

aobj.trans <- AddMetaData(object = aobj, metadata = predicted.labels)

aobj.trans |> DimPlot(group.by = 'predicted.id')

aobj.trans |> DotPlot(c('prediction.score.TRPV3.lo.KC','prediction.score.TRPV3.hi.KC'),
                      scale = F)

aobj.trans |> VlnPlot('prediction.score.max', group.by = 'predicted.id')

# motif analysis ----------
#BiocManager::install(c('motifmatchr','JASPAR2024',
#'TFBSTools','BSgenome.Hsapiens.UCSC.hg38','chromVAR'))
#library(JASPAR2024)
library(TFBSTools)

aobj@meta.data |>
  as_tibble(rownames = '.cell')

pfm <- getMatrixSet(
  x = '~/append-ssd/work/JASPAR2024.sqlite',
  opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)

hm.skin <- aobj |>
  AddMotifs(pfm = pfm,
            genome = BSgenome.Hsapiens.UCSC.hg38)

# cost ~5 min
da_peaks <- aobj |> FindMarkers(
  ident.1 = 10,
  only.pos = TRUE,
  test.use = 'LR',
  min.pct = 0.05,
  latent.vars = 'nCount_peaks',
  assay = 'peaks'
)

da_peaks.neg <- aobj |> FindMarkers(
  ident.1 = 10,
  only.pos = TRUE,
  test.use = 'LR',
  min.pct = 0.05,
  latent.vars = 'nCount_peaks',
  assay = 'peaks'
)

# get top differentially accessible peaks
top.da.peak <- da_peaks |>
  as_tibble(rownames = 'peak') |>
  dplyr::filter(p_val < .005, pct.1 > .2) |>
  pull(peak)

top.da.peak.neg <- da_peaks.neg |>
  as_tibble(rownames = 'peak') |>
  dplyr::filter(p_val_adj < .05, pct.1 > .2) |>
  pull(peak)

# find peaks open in KC
open.peaks <- aobj |>
  AccessiblePeaks(idents = c(2,10))

# match the overall GC content in the peak set
meta.feature <- aobj |>
  GetAssayData(assay = "peaks", layer = "meta.features")

peaks.matched <- MatchRegionStats(
  meta.feature = meta.feature[open.peaks, ],
  query.feature = meta.feature[top.da.peak, ],
  n = 50000
)

# test enrichment
enriched.motifs <- aobj |>
  FindMotifs(
    features = top.da.peak,
    assay = 'peaks',
    #background = peaks.matched
  )

enriched.motifs.neg <- aobj |>
  FindMotifs(
    features = top.da.peak.neg,
    assay = 'peaks',
    #background = peaks.matched
  )

enriched.srebp <- enriched.motifs |>
  dplyr::filter(str_detect(motif.name, 'SREB'))

enriched.motifs |>
  write_csv('mission/FPP/psoriasis/v3h.enrich.motif.csv')

enriched.motifs <-
  read_csv('mission/FPP/psoriasis/v3h.enrich.motif.csv')

## rank-pvalue plot ---------
ranked.motifs <- enriched.motifs |>
  mutate(rank = rank(-fold.enrichment)) |>
  relocate(rank)

ranked.high <- ranked.motifs |>
  filter(str_detect(motif.name, 'SREB') | rank < 11, fold.enrichment > 2,
         p.adjust < .05)

ranked.high

ranked.motifs |>
  filter(fold.enrichment > 1) |>
  ggplot(aes(rank, fold.enrichment)) +
  geom_point(size = .3) +
  geom_point(data = ranked.high, color = 'red2', size = 1) +
  geom_text_repel(data = ranked.high, aes(label = motif.name), size = 2) +
  theme_classic() +
  theme_jpub +
  labs(title = 'Enriched TF motifs in V3-hi KC vs V3-lo KC')

publish_pdf('mission/FPP/psoriasis/motif.rank.fc.pdf')

enriched.motifs |>
  dplyr::filter(p.adjust < .05) |>
  slice_max(fold.enrichment, n = 9) |>
  mutate(motif.name = fct_reorder(motif.name, fold.enrichment)) |>
  ggplot(aes(fold.enrichment, motif.name, fill = -log10(p.adjust))) +
  geom_col() +
  scale_fill_gradient(high = 'red', low = 'black') +
  theme_pubr() +
  labs(title = 'Motif enrichment in scATAC of V3-hi KC vs V3lo KC')

### try volcano plot

filter <- dplyr::filter

enriched.motifs |>
  mutate(fold.enrichment = log2(fold.enrichment)) |>
  mutate(gene = motif.name, avg_log2FC = fold.enrichment, p_val_adj = p.adjust) |>
  plot_pub_volc(group1 = 'v3h kc', group2 = 'v3l kc',
                 highlights = ranked.high$motif.name) +
  labs(title = 'Healthy skin ATAC: V3-hi KC vs V3-lo KC')

aobj |> MotifPlot(motifs = 'MA0828.3')

publish_pdf('mission/FPP/psoriasis/srebf2.motif.seq.pdf')

hm.skin |> DimPlot(group.by = 'manual.main')


## Computing motif activities ---------
hm.skin <- hm.skin |> RunChromVAR(genome = BSgenome.Hsapiens.UCSC.hg38)

DefaultAssay(hm.skin) <- 'chromvar'

# look at the activity of SREBF2
aobj |> FeaturePlot(
  features = enriched.srebp$motif,
  min.cutoff = 'q10',
  max.cutoff = 'q90'
)

## DEGA on chromvar -------
differential.activity <- aobj |>
  FindMarkers(
  ident.1 = 'V3-hi cell',
  group.by = 'manual.main',
  only.pos = TRUE,
  mean.fxn = rowMeans,
  fc.name = "avg_diff",
  assay = 'chromvar'
)

differential.activity |>
  as_tibble(rownames = 'foo')

motif.meta <- aobj |>
  GetMotifData(slot = 'motif.names') 

motif.srebf <- motif.meta |>
  as_tibble() |>
  pivot_longer(everything()) |>
  dplyr::filter(str_detect(value, 'SREBF'))

motif.srebf

MotifPlot(
  object = aobj,
  motifs = 'MA0828.1',
  assay = 'peaks'
)

aobj |>
  DimPlot(group.by = 'manual.main')

aobj |>
  DotPlot(assay = 'chromvar',
          features = motif.srebf$name,
          group.by = 'manual.main',
          cols = 'RdYlBu') +
  RotatedAxis() +
  labs(title = 'SREBF1/2 binding motif openness', x = 'Motif name',
       y = 'Cell type')

aobj |>
  mutate(manual.main = if_else(str_detect(manual.main, 'V3'), 'KC', manual.main)) |>
  DotPlot(assay = 'chromvar',
          features = motif.srebf$name,
          group.by = 'manual.main',
          cols = 'RdYlBu') +
  RotatedAxis() +
  labs(title = 'SREBF1/2 binding motif openness', x = 'Motif name',
       y = 'Cell type')

# motif footprint ----------
aobj.kc <- aobj |>
  filter(str_detect(manual.main, 'V3'))

# take 2 min
# take 1m20s if in.peaks = T
# take 1m40s if change flank from 250bp to 60bp
system.time(mtf.fp <- aobj.kc |>
  Footprint(motif.name = c("SREBF2.1", "ELK1::SREBF2"),
            upstream = 60, downstream = 60,
            genome = BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38)
)

mtf.fp |>
  write_rds('mission/FPP/psoriasis/motif.footprint.rds')

mtf.fp <-
  read_rds('mission/FPP/psoriasis/motif.footprint.rds')

g2 <- mtf.fp |>
  PlotFootprint(c("SREBF2.1", "ELK1::SREBF2"), group.by = 'manual.main',
                normalization = 'subtract') +
  patchwork::plot_layout(ncol = 1)

g2[[1]][[1]] +
  ylim(c(-3,3))

g2[[1]][[1]][["data"]] |>
  filter(flanks) |>
  summarise(max(norm.value), min(norm.value))

g2[[2]][[1]][["data"]] |>
  filter(flanks) |>
  summarise(max(norm.value), min(norm.value))


gg1 <- g2[[1]][[1]][["data"]] |>
  filter(norm.value > -.5) |>
  ggplot(aes(position, norm.value, color = group)) +
  geom_path(linewidth = .2) +
  theme_jpub +
  labs(title = 'SREBF2.1', x = 'Distance from motif',
       y = 'Tn5 insertion\nenrichment')

gg1

gg2 <- g2[[2]][[1]][["data"]] |>
  filter(norm.value > -2) |>
  ggplot(aes(position, norm.value, color = group)) +
  geom_path(linewidth = .2) +
  theme_jpub +
  labs(title = 'ELK1::SREBF2', x = 'Distance from motif',
       y = 'Tn5 insertion\nenrichment')

gg1 / gg2

publish_pdf('mission/FPP/psoriasis/srebf2.footprint.pdf', width = 60)

# coverage plot----------
aobj |>
  CoveragePlot(region = c("TRPV3"),
               group.by = 'manual.main') +
  theme_jpub

aobj |>
  LinkPlot(region = "chr17-3400000-3560000")
